Tunable heat pump by modulating the coupling to the leads 
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We follow the nonequilibrium Green's function formalism to study time-dependent thermal trans- 
port in a linear chain system consisting of two semi-infinite leads connected together by a coupling 
that is harmonically modulated in time. The modulation is driven by an external agent that can 
absorb and emit energy. We determine the energy current flowing out of the leads exactly by solv- 
ing numerically the Dyson equation for the contour-ordered Green's function. The amplitude of 
the modulated coupling is of the same order as the interparticle coupling within each lead. When 
the leads have the same temperature, our numerical results show that modulating the coupling 
between the leads may direct energy to either flow into the leads simultaneously or flow out of the 
leads simultaneously, depending on the values of the driving frequency and temperature. A special 
combination of values of the driving frequency and temperature exists wherein no net energy flows 
into or out of the leads, even for long times. When one of the leads is warmer than the other, 
net energy flows out of the warmer lead. For the cooler lead, however, the direction of the energy 
current flow depends on the values of the driving frequency and temperature. In addition, we find 
transient effects to become more pronounced for higher values of the driving frequency. 

PACS numbers: 05.70.Ln,44.10.+i,63.22.-m,66.70.Lm 



I. INTRODUCTION 



The transport of phonons in mesoscale and nanoscale 
devices is an important issue relevant to the questions 
of heat generation in devices and their structural stabil- 
ity. Experiments measuring the heat generated in elec- 
tric current-carrying metal-molecule junctions found that 
the generated heat can be substantial [l| and can there- 
fore threaten the device's integrity. Efficiently dissipat- 
ing heat in such devices is thus important and a problem 
that must be considered. A way of manipulating heat in 
nanoscale devices is by utilizing a heat pump that directs 
heat from one part of the device to another or to an exter- 
nal reservoir by means of an applied external work. Mod- 
els on the mechanism of such a nanoscale heat pump have 
been proposed in systems where the pump works against 
the thermal gradient between two reservoirs in the sys- 
tem [2|tJ] and in systems where there is no net thermal 
bias between the two reservoirs 0-0|- In addition, other 
models employing quantum particle pumps that differen- 
tiate and filter hot and cold particles have been proposed 
[8j]. In this paper we present an alternative model of a 
phonon pump that directs energy, and thus heat, into 
or out of regions of the device by dynamically modu- 
lating the coupling between those regions. The model 
is different from previous models of heat pumps where 
requirements of either modulating the temperatures of 
reservoirs [g, |7j , or having an external driving force act- 
ing at the central portion of the device [3H5J1, or filtering 
particles according to their temperatures 2, 8] should be 
satisfied. Our model, in comparison, requires an external 



agent that can either absorb or release heat and is dy- 
namically modulating the coupling between two parts of 
the device. A thermal gradient is not necessary for our 
model phonon pump to work. 

To induce a phonon pump action in our model the 
coupling between the two parts of the system is harmon- 
ically modulated in time. Experimentally, this can be 
done by, for example, harmonically varying the distance 
between two molecules, therefore modulating the cou- 
pling between the molecules. Time-dependent transport 
of phonons in molecular systems, however, is a topic that 
is not yet fully understood theoretically. Although our 
understanding of the subject has improved tremendously 
during the past few years, most of the results pertain to 
steady-state and long-time behavior [9j . Time-dependent 
phonon transport with non-adiabatic and strong pertur- 
bations has recently been studied in a thermal switch 
device where the coupling to the reservoirs is abruptly 
turned on [lOj. In this paper we extend the thermal 
switch model to one where the reservoir coupling is 
modulated in time. The system subsequently acts as a 
phonon pump that can be tuned by varying the frequency 
of modulation of the coupling to the reservoirs. 



II. MODEL AND THEORETICAL APPROACH 

In this paper we consider phonon transport in a one- 
dimensional chain. Shown in Fig. [T] is a linear chain con- 
sisting of two semi-infinite leads, or reservoirs, coupled 
together by a coupling k(t) that is modulated in time. 
The particles in each lead are coupled to their nearest 



modulated coupling 

I 




semi-infinite left lead 



semi-infinite right lead 



FIG. 1: (Color online) An illustration of an infinite linear 
chain whose two semi-infinite parts are connected by a cou- 
pling k(t) that is modulated in time. The labels of the first 
three particles in each lead are shown beside each particle. 
Nearest-neighbor particles within the leads interact via an in- 
terparticle spring constant k. An on-site spring with spring 
constant fco is also experienced by each particle. 



neighbors by a coupling constant fc. In addition, each 
particle experiences an on-site potential with spring con- 
stant fen. The temperatures in the left and right leads are 
Tl and Tr, respectively. The particles can vibrate only 
along the horizontal axis. Particles in each lead follow 
a Hamiltonian with purely harmonic, nearest-neighbor, 
interactions that do not vary in time: 

ffQ = jE("") 2 + ^E<^K' « = l,r. (i) 



where the first sum is over all particles and the second 
sum is over all nearest-neighbor pairs in the leads. Each 
nearest-neighbor pair, however, is considered twice and 
so we divide the second term by two. The transformed 
coordinates Ui = \/mxi is used, where x% is the relative 
displacement of the i-th particle of mass m, and K Q is 
the coupling matrix. This matrix is from the dynamic 
matrix of the full system: 



K = 



K L yLR 

yRL K R 



(2) 



where in the one-dimensional chain the spring constant 
matrices K and K are semi-infinite tridiagonal sub- 
matrices consisting of 2fc + fco along the diagonal and — fc 
along both off-diagonals. The coupling matrices V LR (£) 
and V RL (£) are the couplings to the leads that vary in 
time. The Hamiltonian involving the lead coupling is: 



H Lti (t) 






li^(<)«f, 



(3) 



where the sum is over all particles in each lead that are 
directly coupled to the other lead. For the linear chain 
shown in Fig. [1] the coupling matrices each have only one 
nonzero clement: 



C(*) = ~ k (t) and V^{t) = -k(t) 



(4) 



Notice that during the same instant in time, H hR (t) — 
H Rh (t). The total time-dependent Hamiltonian for the 
two-lead system is 



H(t) = H L + H K + H Lti (t) 



(5) 



Note that what we have is an open system consisting of 
the two leads and their coupling that is being modulated 
by an external agent. The energy, therefore, that this 
two-lead system gains or loses is coming from or going to 
the external agent. 

The energy current flowing out of the left lead is 



J L (i) 



dH L 
~df 



([H\H]), 



(6) 



i.e., it is the negative of the expectation value of the rate 
of change in H L . The Heisenberg equation of motion is 
used in the second equality. The position and momentum 
of a particle obey the commutation relation 



u?(t),u?(t)]=ih5 ij 5 al3 , a,/3 = L,R, (7) 



where i and j are particle labels. Note that this commu- 
tation relation only exists at the same instant in time. 
We thus find that the only term in H that does not com- 
mute with H L is H LR . Equation © therefore becomes 



/ L (i) 






L,X ..Lt/LR^R]^ 



mn ^n \ 



(8) 



where all the terms in the right-hand side occur at the 
same time t. Now define the real-time lesser Green's 
function as 



G 



RL,<, 

ij 



(ti,t 2 ) = -Utf(t 2 )uf(t 1 )) 



(9) 



We would like to use this Green's function to properly 
calculate the current. Re-expressing Eq. ([S]) using two 
time variables we get 



jmn 

^(ti)^ R (ti) w R (ii)])| ti=t2=tI (10) 



where we set t\ = t 2 = t in the end. Since V LR (£) = 
V RL (£), we have 



I L (t) 



E C jm(t2,tl) Oil) j- {(u^(t 2 )u^(h)) 



2/i 



j run 



dt 2 



+ Kiti)u){t 2 ))}\ ti=t2=t , 



(ii) 



where Cj m (t2,ti) = [w L (t 2 ),u^(ti)] . Making use of the 
lesser Green's function, Eq. ©, and its complex conju- 
gate, we get 



/ L (i) 



E C ^(*2.*l)^mn(*l) 



jmn 



i m <;A G R L,< (il , i2 ) 



(12) 



where "Im" means taking the imaginary part. Equa- 
tion (1121) is a general equation for the energy current 



flowing out of the left lead. Note that the result is in- 
dependent of the dimension of the system. For the one- 
dimensional linear chain setup shown in Fig. [1] we use 
Eq. (Ql and the fact that Cj m (t,t) = —ihdj m to get 



I L (t)=hk(t)lm\—G w 







RL,< 



(ti,ta) 



(13) 



Equation (p~3|) is the primary equation we use to calculate 
the current flowing out of the left lead. Similarly, we 
define the energy current flowing out of the right lead as 
I R (i) = — (dH K /dt) and derive an expression in terms 
of the lesser Green's function: 



decoupled and each is in thermal equilibrium with tem- 
perature T L and Tr, respectively. At time t = the 
coupling between the leads, k(t), is switched on. We 
then calculate the energy current at time t, i.e., the time 
at the right edge of the contour. Note that instead of 
adiabatic switch-on, the coupling k(t) is abruptly turned 
on at t — 0. In addition, the left and right leads are 
uncorrelated before t = and so there is no imaginary 
tail when the contour goes back to time t = 0. We now 
define the contour-ordered Green's function 



G^(r 1 ,T 2 ) = --(T c uf(r 1 )^(r 2 )) 



(15) 



I R (t) =hk(t)lm 



9 WLR,</, 



fa) 



(14) 



To calculate the currents in Eqs. (fT3j) and (fi"4"|) we 
need to determine the nonequilibrium lesser Green's func- 
tions in those equations. We now follow the Schwinger- 
Keldysh formalism, in which a complex-time contour is 
employed [lll - H^ . to determine the Green's functions. 



where the T c is the contour-ordering operator, T\ and 
r 2 are contour variables, and u R and vh are operators 
in the Heisenberg picture. Transforming to the interac- 
tion picture, we separate the Hamiltonian in Eq. ([5]) into 
the free-particle quadratic part, Hq = H L + H R , and 
the interaction part, H{ n t(t) = H w (t). We then write 
the contour-ordered Green's function in the interaction 
picture as 



^2 



FIG. 2: (Color online) The complex-time contour in the 
Keldysh formalism. The path of the contour begins at time 
t — 0, goes to time t, and then goes back to time t — 0. ri 
and T2 are complex-time variables along the contour. 



Gg L (r 1 ,T 2 ) = -i(T c e-^/o^M^ n R (n)n W r2) \ ) 



where c is the contour shown in Fig. [5] and the average is 
now taken with respect to the equilibrium distributions 
when t < 0. We expand the exponential to perform a 
perturbation expansion. For the Oth-order term we hnd 



G^ (n,r 2 ) = -- (T c uf(n)^(r 2 )) = (17) 



Shown in Fig. [5] is the Keldysh contour we use. At 
times t < we consider the left and right leads to be 



because there is no coupling term that would connect the 
left and right particles. The lst-order term is 



J 



G ij,i(n,T2) 



= (4) T, dT ' V ^(r')(T c ul(T')u K (T')uf(n)^(T 2 )), 

^ ' rnn c 

= E l dT ' f-J(Tc«^(T')«5'(7 5 )>}^(/){-i<T c ««(r')« < R (n)> 



(18) 



where Wick's theorem is used to expand the four-particle 
average into two two-particle averages. There are actu- 
ally three different ways to expand the four-particle av- 
erage but the other two configurations vanish and we are 
left only with the term shown in the second equality. We 
would like to note that the use of Wick's theorem here 
is justified because the expansion is with respect to the 
quadratic Hq. 



Define the equilibrium Green's function of the free 



leads as 



<?g(ri,T 2 ) - -- (T c <(nK(r 2 )) , a = L,R, (19) 



where the average is taken with respect to H , i.e., the 
corresponding equilibrium distribution of the leads. Note 
that unlike the nonequilibrium G RL , the equilibrium g a 
satisfies time-translation invariance and thus its Fourier 
transform exists and can be calculated. Writing Eq. d 



in terms of the equilibrium Green's functions we get 

G$i(ri,7- a ) = X; f^(7i,T')C(r')4(T' I 7i). 

(20) 
To get the lesser version of the nonequilibrium Gfh we 
employ analytic continuation and Langreth's theorem 

14] to get 

mn ^° 

X {^ r (tl,« / )^i < (*'.*2) 

+ ^ < (*i.* , )^ i °(* , ,*2)}, (21) 

where <? in ' r and g m '? are the retarded and advanced ver- 
sions of the equilibrium Green's functions, respectively. 
Similarly, the retarded and advanced versions of the first- 
order nonequilibrium Green's functions are 

G £i' C = £ f d* g%Hti,t')V%;{t')g%(t',t 2 ), (22) 

mn ^° 

where £ = r,a. For the one-dimensional chain shown in 
Fig. [1] only the i = 1 and j = label combination is 
nonzero. In addition, the coupling potential in Eq. (j4]) is 
nonzero only when n = 1 and m = in the sum. All the 
other combinations of the indices do not contribute. 

R a ij R L "0 L R <>' L 



^2 
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X 2 Xi X 2 
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FIG. 3: Diagram representations for the equilibrium Green's 
functions (a) <?5( Tl ' r2 ) an d (b) <?ij(i~i,'>~2), and (c) the 
nonequilibrium Green's function Gy L (n,T2). 

To facilitate the calculation of higher-order terms in 
the perturbation expansion we utilize a diagrammatic 
approach. Shown in Fig. [3] are the diagrams represent- 
ing the relevant contour-ordered Green's functions in the 
expansion. A diagram can not be constructed for the 
zeroth-order term, Eq. (|T7|) . since it contains a u L and 
a u R pair. The first-order term, in contrast, has just 
the right number of double pairs of u L and u R and a 
y RL to connect the two equilibrium Green's functions. 
The second-order term has the same shortcoming as the 
zeroth-order term, i.e., there is an extra it L and u R pair. 
In fact, all the rest of the even-ordered terms have the 
same extra u L and u R pair. All of the even-ordered terms 
therefore do not contribute to the perturbation expan- 
sion. As for the odd-ordered terms, they consist of repe- 
titions of the first-order diagram. Shown in Fig. 3] is the 
perturbation expansion in diagram representation. 

In the first equality of Fig. |4l the first and third- 
order terms are explicitly shown. Notice that the third- 
order term consists of two first-order terms. The dia- 
gram representation of the Dyson equation is shown in 



R R L L 

i , * +j 

Tl T' T 2 

R R L L R R L 
+ '• * X X— 

x, x' x" x'" 

+ 0(v 5 ) 

R R L L 

it * *j 

X, x' T 2 

R R L L R L 

+ t* X X •) 

tl x' t" t 2 



FIG. 4: (Color online) Diagrammatic perturbation expansion 
of G RL in terms of g R and g L . Each interaction vertex, i.e., the 
crossed-out (red) dot, includes an interaction potential V nh 
and an integration with respect to the internal complex-time 
variable. The second equality is the diagram representation 
of the Dyson equation. 



the second equality of Fig. 2J Note that this equation 
includes all terms in the expansion. Based on the dia- 
grammatic rules, we can now write the Dyson equation 
for the nonequilibrium Green's function as 

Gff(Ti,T 2 ) = J2j dT '9f m (n,T')V^(r')glt j (r',r 2 ) 

mn c 

+ E / dT ' J dT " 9fm{rur')V^{r')gl p {r'y') 

mnpq c c 

x v£V')eW,75). (23) 

Applying Langreth's theorem [1J| to Eq. (|23p and then 
iterating the resulting equation, we get a closed-form for- 
mula that includes all orders of the expansion: 



Gi-i ' (h,t 2 ) - G l% { iti,t 2 ) 



£ f dt' Gf^{t x ,t') *£«(f ) G%ftf, t 2 ) 

mn ^° 

E f dt ' GKi^i. *') V^(f) G*f- a (t>, h) 

mn ^° 

^ ft dt> f t dt"Gf^ r (t 1 ,t')V^(t') 



x Gll : <{t',t")V^{t")G^ a {t"M). (24) 

Equation l|2"4"|) is the exact formula for the general 
nonequilibrium Green's function needed to calculate the 
current flowing through the left lead in Eq. flj"2"1) . For 
the linear chain shown in Fig. [1] we use the coupling po- 
tential in Eq. Q and the indices i = n = q = 1 and 
j = m = p = 0. 

To solve Eq. flU} we need to determine Gf^, Gf^' a , 

and G i • ' r . From Eq. (|2ip the first-order nonequilibrium 
Green's function can be calculated by 



G 10 \ (h,t 2 ) - - 



J* dt'k{t') {sS r (*i,05o6 < (*',*a) 



+ Sn' < (^')3o6 a (^2)}- (25) 

To determine the full retarded Green's function we apply 
Langreth's theorem to Eq. ([23)) to get the equation 

G^- r (h,t 2 ) + fdt'k{t')G^{tut')G^' T {t',t2) 
Jo 

(26) 



^RL,r 



^lO.f^l)^), 



where Cl = (w + ir\) 2 — 2k — ko and the choice between 
the plus or minus sign depends on satisfying |A| < 1. 
Furthermore, the lesser equilibrium Green's function can 
be determined from 



15] 



3"' <r 



[u>\=2if a lm{g% r [u]}, _ (29) 

where f a is the Bose-Einstein distribution function of 
the a lead. Given a function F[u>] in frequency space its 
inverse Fourier transform is 



where Eq. (|22|) is used for the first-order Green's func- 
tions. A similar equation can be derived for the full ad- 
vanced Green's function. These two equations for the full 
retarded and advanced Green's functions can be solved 
by the process discussed in Sec. IIIII 



III. NUMERICALLY CALCULATING THE 
ENERGY CURRENT 

In the linear chain, the energy flowing out of the left 
lead, I h (t), at time t can be calculated using Eq. (flUl) . 
the nonequilibrum lesser Green's function G KL,< (ti,t2) 
shown in Eq. (|24[) . and its derivative with respect to 
ti. From Eq. (f24|) . the G RL,< can be calculated from 
the first-order nonequilibrium lesser Green's function 
G 1 ' , the full nonequilibrium retarded and advanced 
Green's functions, G RL ' r and G RL,a , respectively, and 
their derivatives with respect to t 2 . Furthermore, from 
Eq. (|21[) , the G-l ' can be calculated from the integral of 
equilibrium Green's functions g R ' r , ff L,< , g R< , and g h ' a . 
In addition, from Eq. (|26p . the full nonequilibrium re- 
tarded and advanced Green's functions can be calculated 
from the integral of equilibrium Green's functions. All of 
the nonequilibrium Green's functions, therefore, can be 
calculated from the integrals of equilibrium Green's func- 
tions. What we ultimately need then are the equilibrium 
Green's functions. 

Equilibrium Green's functions satisfy time-translation 
invariance and therefore their Fourier transforms exist. 
In frequency space the retarded equilibrium Green's func- 
tions for the semi-infinite linear chain leads are known to 
be M 



9ij 



[ u ] = --\\*-i\, a = L,R, (27) 



where the square brackets mean that the function is a 
Fourier transform and 



^-2^2^^' 



(28) 



F{t!,t 2 ) 



did 

2^ 



F[u\e 



— iuj(t 1 —t2) 



(30) 



We can thus take the inverse Fourier transform of 
Eqs. (|27l) and (|29p to determine the time-dependence of 
the corresponding equilibrium Green's functions. In ad- 
dition, the advanced equilibrium Green's functions can 
be calculated from the retarded version by [9j 



a, a r i 

9a M 



(%'>])* 



(31) 



The integrals appearing in the inverse Fourier transforms 

are numerically calculated using the trapezoidal rule [16| ■ 

After numerically calculating the equilibrium Green's 

functions, we can use Eq. (|25|) to determine the first- 



RL,< 



order nonequilibrium lesser Green's function G i ■ ■[ . We 
again use the trapezoidal rule to calculate the integral in 
Eq. (USD. 

To calculate the full nonequilibrium retarded Green's 
function G i ■ ,r , we solve Eq. (|26|) . This equation is in 
the form of a Frcdholm equation of the second kind [la] 



f{t a ,t b )+ dt' fi(t a ,t')k(t')f(t',t b )=f 1 (t a ,t b ), (32) 



where k and f\ are assumed known and / is the un- 
known. To solve for / we discretize the time into N total 
intervals of incremental length At = t/N. Applying the 
trapezoidal rule to the integral in Eq. (|26|) we get 

N-l 
+ E /l(*a,*j)fc(tj)/(*i,*6) 

+ \h{ta,t N )k{t N )f{t Nl t b )X =fl(t a ,t b ), (33) 

for a set of values of t a and £&. We can recast the calcu- 
lation into a linear problem of the form 



/ f(t ,h) \ 

f(h,t b ) 

\f(tN,t b )J 



At 



( |/i(to,*o)*(*o) 
|/i(ti,to)A(*o) 



fl{tl,tl) 



fc(ti) 



\5/i(<iV,*o)A:(<Jv) /i(tjv,*i)*(*i 



|/i(*o,*jv)*(*jv) \ 

fi(ti,t N )k(t N ) 



? 



/(*i,*b) 



/ /l(t ,*6) \ 
/l(*l,*6) 



\fl{t N ,t N )k{t N ) ) \f(t N ,t b )J \ &&*,%) J 

(Ml 



We end up with a linear problem of the form 

(1 + F) ./ = /i- 



(35) 



The unknown vector / can be determined by decompos- 
ing (1 + F) using LU decomposition and then back sub- 
stituting the result to the vector f\. The G w ' r (ti,t2) 
in Eq. (|26p is numerically determined this way for values 
of ii and ti within the interval [0,t]. The same calcula- 
tion can also be done to determine the advanced Green's 



function G 10 '"(ii,^)- 



Having numerically calculated G 



RL,< 
10,1 ' 



G w ' r , and 



G 



RL,a 
10 ' 
RL,< 



we can then use Eq. (j24|) to determine 



G w ' (ti,<2) and its derivative with respect to t%. The 



current I L (t) is then calculated from Eq. (fl3]) . The same 
steps can be followed to calculate the energy current flow- 
ing out of the right lead, / R (t). 



IV. NUMERICAL RESULTS 

We numerically calculate the time-dependent behavior 
of the energy current flowing out of the leads. The cou- 
pling between the leads is harmonically modulated in the 
form 



k 
k(t) =2^ -cosw d £), 



(36) 



where u>d is the driving frequency. Note that this mod- 
ulated coupling has the same order as k and thus, a 
perturbative calculation is not expected to produce ac- 
curate results. In contrast to perturbative calculations, 
we calculate the current exactly by numerically solving 
the Dyson equation. In all of our calculations, we set 
the interparticle coupling k — 0.625 eV/(A 2 u) and the 
on-site spring constant ko = 0.0625 eV/(A 2 u). These 
choices lead to a natural time scale that we also use as 
the unit of time, [t] = 10~ 14 s. We use a time increment 
of At = 0.1 [t] in our calculations. In addition, choosing 
the values of k and ko a l so se ts the width of the phonon 
band. In a linear chain, the phonon density of states is 
confined to be within ko < lo 2 < 4k + ko. 

We explore several variations of our setup. First, we 
study the energy current when there is no thermal bias 
between the leads. Let the temperature of the leads be 
T = T L = T R . Shown in Fig. [5] is the time-dependent 
behavior of the current flowing out of the left lead for 
four different driving frequencies and three different lead 




FIG. 5: (Color online) The current I L (t) as a function of 
time when the driving frequency is (a) uid = 0.125 [1/t], (b) 
cj d = 0.25 [1/t], (c) uj d = 0.5 [1/t], and (d) uj d = 1 [1/t]. The 
leads have the same temperature T — TL = Tr with values 
T = 10 K (green squares), T — 300 K (red triangles), and 
T = 500 K (blue circles). The harmonic modulation of the 
coupling is shown (red dashed line) as a guide. Its amplitude 
is not drawn to scale. 



temperatures. Notice that the current does not exactly 
follow the modulation. 

Although the amplitude of the modulated coupling is 
kept constant at k/2, as the driving frequency is increased 
the peaks in the current also increases. In addition, as 
the driving frequency becomes sufficiently high, transient 
behavior in the current becomes visible. Transient be- 
havior becomes pronounced when the modulation of the 
coupling changes rapidly. We look at transient behavior 
more closely in Fig. [TO] 

Shown in Fig. [6] are plots of the Fourier transforms of 
the left current when the driving frequencies are Ud — 
0.25 [1/t] and u>d = 1 [1/i]- The peaks occur at frequen- 
cies that are integer multiples of u)d- Modulating the 
lead coupling therefore produces a dynamic energy cur- 
rent that is composed of the first few harmonics of the 
driving frequency. 

In Fig.[5l we notice that the current appears to be more 
on the negative axis as the driving frequency is increased. 
Since we define the left current in Eq. ([6]) as the energy 
flowing out of the left lead, a negative value means that 
the energy is flowing into the lead. To be more definite, 
we calculate how much energy has flowed into the left 
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FIG. 6: (Color online) The square of the Fourier transform 
of the current |/ L [u;]| 2 as a function of the frequency ui when 
the driving frequencies are (a) uid = 0.25 [1/t] and (b) iQd = 
1 [1/t]. The leads have the same temperature T and shown 
are when T = 100 K (A, red) and T = 300 K • blue). 
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FIG. 7: (Color online) The energy E h (t) as a function of 
time for driving frequencies (a) ujd = 0.125 [1/t], (b) w<j = 
0.25 [1/t], (c) uj c1 = 0.5 [1/t], and (d) uj d = 1 [1/t]. The leads 
have the same temperature T and shown are for T = 10 K 
(green squares), T — 300 K (red triangles), and T = 500 
K (blue circles). The harmonic modulation of the coupling 
is shown (red dashed line) as a guide. Its amplitude is not 
drawn to scale. 



between the leads, the left and right leads are indistin- 
guishable. The plots for i R (i) and E R (t) are therefore 
exactly the same as those shown for the left lead in Fig. [5] 
and Fig. [JJ respectively. The external agent therefore 
supplies energy to both the left and right leads. 

Decreasing the driving frequency, we see from Fig. [JJ 
that for certain values of the temperature, instead of the 
energy flowing into the leads, energy is actually flowing 
out from the leads. In Fig. [7[c), for example, when the 
driving frequency is cud = 0.5 [1/t], energy flows into the 
leads when T = 300 K and T = 500 K but it flows out of 
the leads when T = 10 K. 
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FIG. 8: (Color online) The energy E L as a function of the 
lead temperature T, which is the same for both leads. Each 
curve in the plot corresponds to a specific value of time, (a) 
For lj = 0.25 [1/t] at time t = 11.8 [t] (T, red), t = 36.9 [t] 
(■, green), t = 62.1 [t] (A, orange), and t = 87.2 [t] (•, blue), 
(b) For oj = 0.5 [1/t] at time t = 19.4 [t] (T, red), t = 44.6 [t] 
■ green), t = 69.7 [t] (A, orange), and t = 94.8 [t] (•, blue). 



lead by 



E L (t) = / I L (t')dt'. 



(37) 



We again use the trapezoidal rule to numerically calcu- 
late the integral. Shown in Fig. [7] is the energy E h (t) 
for the four different driving frequencies and three differ- 
ent temperatures corresponding to those shown in Fig. [5] 
Notice that for the highest driving frequency, as shown 
in Fig. [TJd), the energy increases negatively in time for 
all three temperatures shown. Therefore, harmonically 
modulating the lead coupling by a driving frequency of 
ujd = 1 [1/t] moves energy from the external agent into 
the left lead. Furthermore, since there is no thermal bias 



Shown in Fig. [S] are plots of how E L varies for different 
values of the driving frequency and lead temperature, at 
specific instants of time. Note that in Fig. [JJthe energy 
E L oscillates in time. The data points in Fig.[8]are chosen 
from Fig. [7J during the times when E L is at a trough 
in the oscillating energy. We can also choose different 
sets of data points, e.g., points at the crest instead of 
the trough, and then plot them like those in Fig. [8] The 
plots, however, will look the same except for a translation 
along the vertical axis. 

In Fig. [8l the E L values are negative at higher tem- 
peratures. At lower temperatures, however, E h can be 
positive, depending on the values of ujd and t. Notice 
that there is a temperature T c where the curves, for 
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one value of uid, intersect. In Fig. [8ja), for example, 
T c rj 300 K. This T c is lower for higher lo^ values, as 
shown in Fig. |UJb). Now, since the choice of taking data 
points only at the trough of E h is arbitrary, we may also 
choose a different location along E L such that at T c we 
get E^(T C ) = 0, for any time t. For such a choice, the 
values of E h below T c are positive and above T c are neg- 
ative. 

Notice in Fig. [8] that the slope of the curves becomes 
steeper at later times. The value of E L at T = T c , how- 
ever, remains the same at any time t. Therefore, at much 
later times when the slope of the curves are much steeper, 
when T > T c both the left and right leads absorb energy 
from the external agent, while when T < T c the leads 
emit energy to the external agent. 
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FIG. 9: (Color online) Time evolution of (a) the current in the 
left lead, I L , (b) the right lead, I R , (c) the sum of the currents, 
I B = I -\-I , and (d) the sum of the energies in both leads, E s . 
The average temperature between the leads are T avc = 10 K 
(green squares), T aV G = 300 K (red triangles), and T avc = 500 
K (blue circles). The driving frequency uid = 0.25 [1/t]. The 
harmonic modulation of the coupling is shown (red dashed 
line) as a guide. Its amplitude is not drawn to scale. 



We now consider what happens when the leads have 
different temperatures. Let the lead temperatures be 



T L = (l + e)T ave , 
2r = (1 — e)T avo , 



(38) 



i.e., the left lead is warmer than the right lead. Heat 
would therefore spontaneously flow, except for transient 
effects, from the left lead to the right lead if the cou- 
pling between the leads is not modulated in time [10j. 
In our calculations we consider a temperature variation 
of £ = 0.1. Figure HJa) shows the current flowing out 
of the left lead as a function of time when the driving 
frequency uid = 0.25 [1/t]. The time evolution of the 
current does not exactly follow the modulation of the 
lead coupling. Taking the Fourier transform of the cur- 
rent produces peaks at a few integer multiples of the 
driving frequency. This situation is similar to the case 
when the leads have the same temperature, as shown in 



Fig. [6] Notice in Fig. HJa) that the current is mostly 
positive and thus, as expected, energy is flowing out of 
the warmer left lead. Shown in Fig. (S^b) is the current 
flowing out of the right lead. For the cooler right lead, 
the question of whether the current is mostly negative or 
positive depends on the temperature of the lead. When 
Tavc = 10 K, the current is mostly positive and therefore, 
energy is mostly flowing out of the cooler right lead. At 
Tkvo = 10 K therefore, the current is flowing out of both 
the left and right leads, i.e., energy from both leads is 
being absorbed by the external agent. In contrast, when 
Tkvc = 500 K the current in the right lead is mostly nega- 
tive and therefore energy is mostly flowing into the cooler 
right lead. 

Notice that the plots in Fig. Eth) are not mirror reflec- 
tions, with respect to the horizontal axis, of the plots in 
Fig. HJa). When the sum of the currents are taken, i.e., 
I s (t) = I L {t) + I K (i), the results are the plots shown in 
Fig. Hfc). Note that although the plots only show the 
time evolution of I s up to time t = 50 [t] our calcula- 
tions extend until time t = 100 [i\. Comparing Fig. [3Jc) 
to Fig. [5jb) , we find the plots to coincide (note that the 
plot in Fig. [5jb) is only for the left lead and so the values 
of the current should be multiplied by two). Further- 
more, we find that the values of I s (t) when e = 0.1 are 
numerically very close to, if not the same as, the values of 
I s (t) when e — at each corresponding instant of time t. 
Our numerical results thus hint that the sum of currents 
I s (£) may be independent of the value of e. 

Modulating the lead coupling when e = 0.1 therefore 
results in plots of the net current I s that numerically 
coincide with those when e = 0. The direction of the 
current flow in each lead, however, differs depending on 
whether the leads have the same or different tempera- 
tures. When the leads have the same temperature the 
current can either flow into both of the leads at the same 
time or flow out of the leads at the same time. In con- 
trast, when the leads have different temperatures, the 
energy can flow out of both the warmer left lead and the 
cooler right lead resulting in a net current that flows out 
of the linear chain system and into the external agent. 
This happens, for example, when T ave = 10 K in Fig. |H1 
Increasing the average temperature to T aV c = 300 K, we 
find a balance between the energy that flows out of the 
warmer left lead and the energy that flows into the cooler 
right lead resulting in no net energy flow for the whole 
system, as shown in Fig. [9jd). Increasing the average 
temperature further to T avo = 500 K, we find that en- 
ergy flows out of the left lead and flows into the cooler 
right lead resulting in a net energy flowing into the chain 
system. We thus see that the current can either flow into 
or out of the cooler right lead depending on the value of 
the average temperature of the leads. 

The driving frequency in Fig. [5] is uid = 0.25 [1/t]. As 
shown in Fig. [8] there is a temperature T c where there 
is no net energy flowing into or out of the linear chain 
system. When uid is varied the T c also varies. Similarly, 
for the case when one of the leads is warmer than the 



other, when T ave < T c we find that the current flows out 
of the cooler right lead resulting in a net energy flowing 
out of the chain system. When T avc > T c we find the 
current to flow into the cooler right lead resulting in a 
net energy flowing into the chain system. Note that for 
any value of T avc energy flows out of the warmer left lead. 



pends on the values of the driving frequency LOd and the 
temperature T of the leads. The faster ix>d and higher T 
produce larger transient current amplitudes. 



SUMMARY 
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FIG. 10: (Color online) The current I L (i) as a function of 
time t when the coupling between the leads is gradually in- 
creased with driving frequencies (a) u)d = 0.125 [1/t], (b) 
u)d = 0.25 [1/t], and (c) uJd — 0.5 [1/t]. Both leads have the 
same temperature T = 10 K (green squares), T = 300 K (red 
triangles), and T = 500 K (blue circles). The gradually in- 
creasing coupling is shown (red dashed line) as a guide. Its 
amplitude is not drawn to scale. 

We now investigate the effects of the speed of modu- 
lation on the transient behavior of the current. In Fig. [5] 
we find that the transient becomes visible as the driving 
frequency u>d is increased. To clearly see the effects of 
how fast the coupling is changing, we consider gradually 
increasing the coupling in the form 



k(t) = k tanhtOdt. 



(39) 



Shown in Fig. [TU] are plots of the energy current in time 
for various values of the driving frequency u>d and temper- 
ature T. We consider the leads to have the same temper- 
ature T and thus, for later times, we expect there to be 
no steady-state current. In Fig. [TJJJwe see that at earlier 
times the transient behavior shows up as rapid bumps in 
the current and then eventually subsides down to zero at 
later times. The amplitude of the transient current de- 



Dynamically modulating the coupling between the 
leads in the form shown in Eq. (|36|) can result in the 
energy current to either flow into or out of the leads, de- 
pending on the values of the driving frequency u>d and 
the lead temperature T, even when the leads have the 
same temperature. For such a case, it is possible for the 
energy current to either flow out of both leads at the 
same time or flow into both leads at the same time, as 
shown in Fig. [7] In addition, in Fig. [S] we see that for a 
given value of u>d there exists a temperature T c where, in 
the long-time limit, when the temperature of the leads 
T < T c we find the current to flow out of both leads 
and when T > T c we find the current to flow into both 
leads. For the case when the leads have different tem- 
peratures, with e = 0.1, the direction of the flow of the 
energy current in the cooler lead depends on the values 
of Tavc and Ud- When T ave < T c the current flows out 
of the cooler lead but when T avo > T c the current flows 
into the cooler lead. Current flows out of the warmer 
lead for any temperature T avo . Gradually increasing the 
lead coupling in the form shown in Eq. (|3T)1) shows that 
the amplitude of the transient depends on how fast the 
lead coupling changes. Faster changes in the lead cou- 
pling result in larger transient amplitudes, as shown in 
Fig. 1101 As a consequence, harmonically modulating the 
lead coupling with a faster driving frequency results in 
a more pronounced transient behavior in the current, as 
shown in Fig. [5] 
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